Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/a_utils/expanded.lha/expanded/src/SubSet.C

C/C++ Source or Header | 1992-03-19 | 44.6 KB | 1,592 lines

//
// Linear-Affine-Projective Geometry Package
//
// SubSet.C
//
// $Header$
//
// William J.R. Longabaugh
// University of Washington
//
// Implementation of the linear-affine-projective geometry
// package described in William J.R. Longabaugh, "An Expanded
// System for Coordinate-Free Geometric Programming", Master's
// thesis, University of Washington, 1992.
//
// Copyright (c) 1992, William J.R. Longabaugh
//   Copying, use and development for non-commercial purposes permitted.
//   All rights for commercial use reserved.
//   This software is unsupported and without warranty; it is
//   provided "as is".
//
// ***********************************************************************

#include <math.h>
#include "Lap.h"

// ***********************************************************************
//
// Class for exclusive use of SubSet class to hold onto data for a
// subset.  What is really wanted is a simple structure; I am not
// really interested in data hiding here.  However, under g++ 1.37.1,
// this apparently MUST be a class, not a struct, so that virtual
// ptrs of class members are initialized.  Also, the fact that class
// members will be automatically initialized using null constructors
// is desirable.  Thus, this is a class, but note that all members
// are public, allowing unlimited access.
// Note that since there is a VSubSet in this block, SubSets
// must in fact be implemented as pointers to information blocks:
//

class SubSetInfo {

  public:

    SubSetType t;                // Type of subset
    char     name[MAX_NAMESIZE];        // Name of subset
    Space    s;                         // Space containing subset
    int      d;                // Dimension of subset
    Matrix   tstmtx;            // Testing matrix
    Matrix   fullbas;            // Matrix rep. of full basis
    int      isOffset;            // Column for testing
    int      substart;            // Column for testing
    int      zerostart;            // Column for testing
    Scalar   offset;            // Offset of affine hyperplane
    Matrix   fromfull;            // Matrix from full to subset basis
    Matrix   tofull;            // Matrix from subset to full basis
    GeObType accepts;            // Type of object in subset
    VSubSet  tansub;            // Tangent subset

// Constructors/Destructors:

        SubSetInfo(void)  {}         // Null constructor
        ~SubSetInfo(void) {}         // Destructor
        SubSetInfo(char* namein, int d); // Used to build error info block
};

// ***********************************************************************
//
// Declare the existence of the error info block:

static SubSetInfo errinfo;


// ***********************************************************************
// ***********************************************************************
//
// SubSetInfo Class
//
// ***********************************************************************
// ***********************************************************************
//
// Special constructor used to build an "error info block".  Uninitialized
// subsets point to that block.

SubSetInfo::SubSetInfo(char* namein, int dim)
{
  t = NULL_SUBSET;
  d = dim;
  isOffset = -1;
  substart = -1;
  zerostart = -1;
  offset = 0.0;
  accepts = NULL_GEOB;

// Local copy of the debug name:

  strncpy(name, namein, MAX_NAMESIZE - 1);
  name[MAX_NAMESIZE - 1] = '\0';
}

// ***********************************************************************
// ***********************************************************************
//
// SubSet Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Great candidates for inlining, but see the comment in the file
// Space.C:
//
// ***********************************************************************

int SubSet::IsOffset(void) {return info->isOffset;}

// ***********************************************************************

Scalar SubSet::Offset(void) {return info->offset;}

// ***********************************************************************

Matrix& SubSet::FromFull(void) {return info->fromfull;}

// ***********************************************************************

Matrix& SubSet::ToFull(void) {return info->tofull;}

// ***********************************************************************

Matrix& SubSet::AugBasis(void) {return info->fullbas;}

// ***********************************************************************
//
// Used to normalize the representation of projective points contained
// in an affine subset.  It does not do ANY error checks, so beware!
//

GeOb SubSet::Normalize(GeOb& ob)
{
  Matrix newmat = ob.MatrixRep();
  Scalar factor = info->offset / (newmat * info->tstmtx)[0][info->isOffset];

  GeOb retval(ob.Holds(), ob.SpaceOf(), factor * newmat);
  return retval;
}

// ***********************************************************************
//
// Given two affine subsets in projective spaces, determine the
// scalar factor that must be applied to place the origin point
// of the first at the offset of the secong origin point
// It does not do ANY error checks, so beware!
//

Scalar SubSet::NormVal(SubSet& targ)
{
  Matrix newmat = targ.ToFull()[targ.info->isOffset];
  Scalar factor = info->offset / (newmat * info->tstmtx)[0][info->isOffset];
  return factor;
}

// ***********************************************************************
//
// This function copies only a maximum number of characters into
// the debug name buffer.
//

void SubSet::CopyDebugName(char* buf)
{
  strncpy(info->name, buf, MAX_NAMESIZE - 1);
  info->name[MAX_NAMESIZE - 1] = '\0';
  return;
}

// ***********************************************************************
//
// This function builds a tagged debug name.  The user is responsible
// for calling DestroyTagName() when finished.
//

char* SubSet::BuildTagName(char* tag, Space& s)
{
  int len = strlen(tag);
  char* buf = new char[MAX_NAMESIZE + len];

  strncpy(buf, tag, len);
  s.Name(buf + len);
  return (buf);
}

// ***********************************************************************
//
// This function builds a tagged debug name.  The user is responsible
// for calling DestroyTagName() when finished.
//

char* SubSet::BuildTagName(char* tag, SubSet& s)
{
  int len = strlen(tag);
  char* buf = new char[MAX_NAMESIZE + len];

  strncpy(buf, tag, len);
  s.Name(buf + len);
  return (buf);
}

// ***********************************************************************
//
// This function kills a tagged debug name.
//

void SubSet::DestroyTagName(char* buf)
{
  delete buf;
}

// ***********************************************************************
//
// This little routine is used by subset constructors to convert lists
// of arbitrary geometric objects into objects of the desired type
// in the desired space.  A matrix containing std. coords of the
// objects is built:
//

Boolean SubSet::ConvertList(GeObList& v, GeObType targ,
                Space& destspace, Matrix* temp)
{
  for (int i = 0; i < v.Length(); i++) {
    GeOb hold = v[i].MapTo(targ);
    if (hold.SpaceOf() != destspace) {
      return (FALSE);
    }
    (*temp)[i] = hold.MatrixRep()[0];
  }
  return (TRUE);
}

// ***********************************************************************
//
// Modified Gram-Schmidt orthogonalization.  The "have" matrix has
// orthonormal rows, and represents the existing orthonormal basis
// of some subspace.  The "try" matrix contains vectors to use
// for constructing a total of "dim" orthonormal vectors; any
// vectors in the try matrix that are not linearly independent
// are skipped.  We set the error flag if there are not enough vectors
// to use for meeting the "dim" target:
//

Matrix SubSet::GSO(Matrix& have, Matrix& try, int start,
           int dim, Boolean* error)
{
  Matrix retval(dim, have.Columns());
  int count = 0;

// Start by filling in with the vectors we already have:

  for (int i = 0; i < have.Rows(); i++) {
    retval[i] = have[i];
    count++;
  }

// Go through the test vectors until we have enough:

  int tryrow = start;
  while (count < dim) {
    Matrix nextrow(1, have.Columns());
    if (tryrow == try.Rows()) {
      *error = TRUE;
      return retval;
    }
    if (this->Orth(&retval, count, try, tryrow)) {
      count++;
    }
    tryrow++;
  }
  *error = FALSE;
  return retval;
}

// ***********************************************************************
//
// Support for GSO.  Does the actual GSO for a vector.  Returns TRUE
// iff the vector could be orthogonalized.  FALSE means it was
// contained in the subspace we have built.
//

Boolean SubSet::Orth(Matrix* have, int count, Matrix& try, int tryrow)
{
  Scalar factor;
  Scalar mag;

// Modified GSO does orthogonalization incrementally, once for each
// vector we already have.  If we ever end up with the zero vector,
// the attempt has failed and we need to return.  Otherwise, normalize
// as we go:

  Matrix hold = try[tryrow];
  for (int i = 0; i < count; i++) {
    Matrix proj = (*have)[i];
    factor = (proj * Transpose(hold))[0][0] / (proj * Transpose(proj))[0][0];
    hold -= (factor * proj);
    mag = sqrt((hold * Transpose(hold))[0][0]);
    if (fabs(mag) < EPSILON) {
      return (FALSE);
    } else {
      hold = hold / mag;
    }
  }
  (*have)[count] = hold[0];
  return (TRUE);
}

// ***********************************************************************
//
// Build a subset from an info block:
//

SubSet::SubSet(SubSetInfo* ts)
{
  info = ts;
  type = ANY_SUBSET;
  holds = ts->t;
}

// ***********************************************************************
//
// Dumps data for a subset:
//

void SubSet::data_out(ostream &c, int indent, char* name)
{
  char *ibloc = new char[indent + 1];
  for (int i = 0; i < indent; i++) {
    *(ibloc + i) = ' ';
  }
  *(ibloc + indent) = '\0';

  c << ibloc << ast;
  c << ibloc << name;
  c << ibloc << "Type of subset: ";
  SubSetTypeOut(c, type);
  c << "\n";
  c << ibloc << "Currently holds: ";
  SubSetTypeOut(c, holds);
  c << "\n";
  if (holds != NULL_SUBSET) {
    c << ibloc << "Name of subset: " << info->name << "\n";
    c << ibloc << "Embedding space:\n";
    info->s.debug_out(c, indent + ERR_IND);
    c << ibloc << "Dimension: " << info->d << "\n";
    c << ibloc << "Type of object in subset: ";
    GeObTypeOut(c, info->accepts);
    c << "\n";
    c << ibloc << "Offset value column: " << info->isOffset << "\n";
    c << ibloc << "Subset value starting column: " << info->substart << "\n";
    c << ibloc << "Zero value starting column: " << info->zerostart << "\n";
    c << ibloc << "Offset of affine hyperplane: " << info->offset << "\n";
    c << ibloc << "Memory location: " << (long)info << "\n";
    c << ibloc << "Matrix representation of testing matrix:\n";
    info->tstmtx.debug_out(c, indent + ERR_IND);
    c << ibloc << "Matrix representation of full basis matrix:\n";
    info->fullbas.debug_out(c, indent + ERR_IND);
    c << ibloc << "Matrix to convert to subset basis from full basis:\n";
    info->fromfull.debug_out(c, indent + ERR_IND);
    c << ibloc << "Matrix to convert to full basis from subset basis:\n";
    info->tofull.debug_out(c, indent + ERR_IND);
    c << ibloc << "Tangent subset:\n";
    info->tansub.debug_out(c, indent + ERR_IND);
  }
  c << ibloc << ast;

  delete ibloc;
  return;
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************

SubSet::SubSet(void) {type = ANY_SUBSET; holds = NULL_SUBSET; info = &errinfo;}

// ***********************************************************************
//
// Memberwise initialization:
//

SubSet::SubSet(SubSet &s)
{
  info = s.info;
  type = ANY_SUBSET;
  holds = s.holds;
}

// ***********************************************************************
//
// Memberwise assignment:
//

SubSet& SubSet::operator=(SubSet &s)
{
//
// This weird checking is required to bypass the apparent inheritance of
// assignment under g++ 1.37:
//
  if ((type != ANY_SUBSET) &&
      ((type != s.holds) && (s.holds != NULL_SUBSET))) {
    errh.ErrorExit("SubSet& SubSet::operator=(SubSet &)",
                   "Illegal assignment attempted", *this, s);
  }
  holds = s.holds;
  info = s.info;
  return (*this);
}

// ***********************************************************************
//
// Output for debugging:
//

void SubSet::debug_out(ostream &c, int indent)
{
  static char* name = "SubSet object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
// Great candidates for inlining, but see the comment in the file
// Space.C:
//
// ***********************************************************************

char* SubSet::Name(char *buf) {return (strcpy(buf, info->name));}

// ***********************************************************************

Space SubSet::EmbeddingSpace(void) {return (info->s);}

// ***********************************************************************

GeObType SubSet::Accepts(void) {return (info->accepts);}

// ***********************************************************************

int SubSet::Dim(void) {return (info->d);}

// ***********************************************************************
//
// Returns TRUE iff the subset spans the entire embedding space.
//

Boolean SubSet::IsFullSpace(void)
{
  if ((info->t == PROJECTIVE_SUBSET) && (info->s.Holds() == VEC_SPACE)) {
    return (info->d + 1 == (info->s).Dim());
  } else if ((info->t == AFFINE_SUBSET) && (info->s.Holds() == PROJ_SPACE)) {
    return (FALSE);
  } else {
    return (info->d == (info->s).Dim());
  }
}

// ***********************************************************************
//
// Returns true iff the subset is projective and was constructed
// by specifying a set of base points.

Boolean SubSet::HasRemovedPoints(void)
{
  return ((holds == PROJECTIVE_SUBSET) && (info->substart != 0));
}

// ***********************************************************************
//
// Returns a list of points that spans the points at infinity.  Only
// applicable for affine subsets defined on projective spaces.

GeObList SubSet::AtInfinity(void)
{
  static char* sig = "GeObList SubSet::AtInfinity(void)";

  if ((holds != AFFINE_SUBSET) || (info->s.Holds() != PROJ_SPACE)) {
    errh.ErrorExit(sig, "Must be affine subset in projective space", *this);
  }

// The tofull matrix encodes the points at infinity.  Extract the
// std. coords. and build the list.

  int rows = info->tofull.Rows() - 1;
  GeObList retval(rows);

  for (int i = 0; i < rows; i++) {
    retval[i] = GeOb(info->accepts, info->s, info->tofull[i]);
  }
  return (retval);
}

// ***********************************************************************
//
// Returns the tangent subset associated with an affine subset.
// If embedding space is affine, the subset is a vector subspace of
// the tangent space.  If embedding space is a vector space, the subset
// is a vector subset of that space.  If it is in a projective
// space, it is an error to use this.
//

VSubSet SubSet::TangentSub(void)
{
  static char* sig = "VSubSet SubSet::TangentSub(void)";

  if (holds != AFFINE_SUBSET) {
    errh.ErrorExit(sig, "Only allowed for affine subsets", *this);
  }
  if ((info->s).Holds() == PROJ_SPACE) {
    errh.ErrorExit(sig, "Embedding space cannot be projective", *this);
  }

// Perhaps the user has already requested the tangent space, so one has
// been built.  Reuse this tangent space.

  if ((info->tansub).Holds() != NULL_SUBSET) {
    return (info->tansub);
  }

// If the embedding space is an affine space, and the subset is the
// same dimension, the tangent subset is simply the full standard
// subset of the tangent space:

  if (((info->s).Holds() == AFF_SPACE) && (this->IsFullSpace())) {
    info->tansub = (info->s).GetSpace(TANGENT).FullSet();
    return (info->tansub);
  }

// Need to build a new tangent space.

  SubSetInfo* temp = new SubSetInfo;
  Boolean isvs = (info->s).Holds() == VEC_SPACE;

  char* tag = this->BuildTagName("Subset tangent to ", *this);
  strncpy(temp->name, tag, MAX_NAMESIZE - 1);
  temp->name[MAX_NAMESIZE - 1] = '\0';
  this->DestroyTagName(tag);

  temp->t = LINEAR_SUBSET;
  temp->s = (isvs) ? info->s : (info->s).GetSpace(TANGENT);
  temp->d = info->d;
  temp->offset = 0.0;
  temp->accepts = (isvs) ? info->accepts : AFF_VECTOR;
  temp->substart = 0;

// If the parent space is a vector space, the situation is simple.
// The tofull matrix omits the offset row, and we expect the
// component in the offset direction to be zero.  Fromfull has
// one less column:

  if (isvs) {
    temp->tofull = Matrix(info->d, info->s.Dim());
    for (int i = 0; i < info->d; i++) {
      temp->tofull[i] = info->tofull[i];
    }
    temp->fullbas = info->fullbas;
    temp->tstmtx = info->tstmtx;
    temp->zerostart = info->isOffset;
    temp->isOffset = -1;
    Matrix trg = ZeroMatrix(temp->s.Dim(), temp->d);
      for (int j = 0; j < temp->d; j++) {
        trg[j][j] = 1.0;
      }
    temp->fromfull = temp->tstmtx * trg;
  } else {

// If the parent space is an affine space, the situation is more involved,
// since we need to build the subset in the separate tangent vector
// space.  Fortunately, the basis for the subset is expressed using
// tangent space vectors and an offset, so all we have to do is
// omit the offset and convert the representation:

    temp->tofull = Matrix(info->d, info->s.Dim());
    Matrix convert = info->tofull * Transpose(info->s.HoistTangent());
    int slot = 0;
    for (int i = 0; i < convert.Rows(); i++) {
      if (i != info->isOffset) {
        temp->tofull[slot++] = convert[i];
      }
    }
    Boolean errflag;
    temp->fullbas = this->GSO(temp->tofull, IdentityMatrix(temp->s.Dim()),
                  0, temp->s.Dim(), &errflag);
    if (errflag) {
      errh.ErrorExit(sig, "Unexpected error", *this);
    }
    temp->tstmtx = Inverse(temp->fullbas);
    Matrix trg = ZeroMatrix(temp->s.Dim(), temp->d);
    for (int k = 0; k < temp->d; k++) {
      trg[k][k] = 1.0;
    }
    temp->fromfull = temp->tstmtx * trg;
    temp->zerostart = info->isOffset;
    temp->isOffset = -1;
  }

// Construct the new subset and return it;

  info->tansub = SubSet(temp);
  return (info->tansub);
}

// ***********************************************************************
//
// Returns true if the given geometric object is in the subset.
// (Actually, it is true if the object is within EPSILON).

Boolean SubSet::IsIn(GeOb &g)
{
  Scalar testval = EPSILON;

// First check if object is in the embedding space.  If not, return
// false:

  GeOb temp = g.MapTo(info->accepts);
  if (temp.SpaceOf() != info->s) {
    return (FALSE);
  }

// Then multiply the object matrix by the testing matrix, and
// compare the coordinate results with the targets.

  Matrix result = temp.MatrixRep() * info->tstmtx;

// Affine subsets need a coord. == offset.  If this is an affine
// subset of a projective space, we just need to make sure the
// offset is not zero (but also scale the epsilon to take into account
// the scaling of the offset to the desired value).  If this is a
// projective subset, we need to make sure the GeOb is not a base point:

  if (info->t == AFFINE_SUBSET) {
    if (info->s.Holds() == PROJ_SPACE) {
      if (fabs(result[0][info->isOffset]) <= EPSILON) {
        return (FALSE);
      } else {
        Scalar factor = fabs(info->offset / result[0][info->isOffset]);
        testval = testval / factor;
      }
    } else {
      if (fabs(result[0][info->isOffset] - info->offset) > EPSILON) {
        return (FALSE);
      }
    }
  } else if (info->t == PROJECTIVE_SUBSET) {
    Boolean allzero = TRUE;
    for (int i = info->substart; i < info->zerostart; i++) {
      if (fabs(result[0][i]) > EPSILON) {
        allzero = FALSE;
    break;
      }
    }
    if (allzero) {
      return (FALSE);
    }
  }

// Coords. for space orthogonal to subset == 0:

  for (int i = info->zerostart; i < result.Columns(); i++) {
    if (fabs(result[0][i]) > testval) {
      return (FALSE);
    }
  }

  return (TRUE);
}

// ***********************************************************************
//
// Given another subset, this function returns TRUE iff the given subset
// is a subset of this subset.  Currently not implemented for
// projective subsets with removed base points.  Note that only
// comparisons of subsets of identical type are allowed.  (The system
// does not test if an affine subset is a subset of a vector subset,
// but simply answers FALSE).

Boolean SubSet::IsSubset(SubSet &s)
{
  static char* sig = "Boolean SubSet::IsSubset(SubSet&)";
  int toprow;

// A quick kill is if the subsets are the same definition:

  if (info == s.info) {
    return (TRUE);
  }

// Check that the operation is valid:

  if (this->HasRemovedPoints() || s.HasRemovedPoints()) {
    errh.ErrorExit(sig,
      "Comparing projective subsets with removed points not implemented",
      *this, s);
  }

// Make sure the embedding spaces and subset types match:

  if (s.info->s != info->s) {
    return (FALSE);
  }
  if (s.holds != holds) {
    return (FALSE);
  }

// Hit the representation with the testing matrix:

  Matrix rep = s.ToFull() * info->tstmtx;

// For affine subsets, the origin of the argument subset frame
// must be in this subset, and the tangent space elements cannot
// have an offset:

  if (s.Holds() == AFFINE_SUBSET) {
    Scalar testval = EPSILON;
    toprow = rep.Rows() - 1;
    if (info->s.Holds() == PROJ_SPACE) {
      if (fabs(rep[s.info->isOffset][info->isOffset]
      * s.info->offset) <= EPSILON) {
        return (FALSE);
      } else {  // Scale test value to reflect normalization to offset
        Scalar factor = fabs(info->offset /
                 rep[s.info->isOffset][info->isOffset]);
        testval = testval / factor;
      }
    } else {
      if (fabs((rep[s.info->isOffset][info->isOffset]
      * s.info->offset) - info->offset) > EPSILON) {
      return (FALSE);
      }
    }
    for (int j = info->zerostart; j < rep.Columns(); j++) {
      if (fabs(rep[s.info->isOffset][j]) > testval) {
        return (FALSE);
      }
    }
    for (int i = 0; i < toprow; i++) {
      if (fabs(rep[i][info->isOffset]) > EPSILON) {
        return (FALSE);
      }
    }
  } else {
    toprow = rep.Rows();
  }

//  Each basis vector for the argument subset or subset tangent
//  space must be contained in this subset or subset tangent space:

  for (int i = 0; i < toprow; i++) {
    for (int j = info->zerostart; j < rep.Columns(); j++) {
      if (fabs(rep[i][j]) > EPSILON) {
        return (FALSE);
      }
    }
  }
  return (TRUE);
}

// ***********************************************************************
// ***********************************************************************
//
// VSubSet Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Builds a standard subset for a vector space.  CAUTION:  Although this
// is a single-argument constructor, it is NOT intended to be used to
// cast a space into a subset!

VSubSet::VSubSet(VSpace& s)
{
// Create the new SubSetInfo:

  info = new SubSetInfo;
  type = LINEAR_SUBSET;
  holds = LINEAR_SUBSET;

  char* tag = this->BuildTagName("Standard subset for ", s);
  this->CopyDebugName(tag);
  this->DestroyTagName(tag);

  info->s = s;
  info->t = holds;
  info->d = s.Dim();
  info->tstmtx = IdentityMatrix(info->d);
  info->fullbas = IdentityMatrix(info->d);
  info->isOffset = -1;
  info->substart = 0;
  info->zerostart = info->d;
  info->offset = 0.0;
  info->fromfull = IdentityMatrix(info->d);
  info->tofull = IdentityMatrix(info->d);
  info->accepts = s.NativeType();
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Memberwise assignment:
//

VSubSet& VSubSet::operator=(VSubSet &s)
{
  info = s.info;
  holds = s.holds;
  return *this;
}

// ***********************************************************************
//
// A subset can only be cast down to a vector subset if it is one (no
// conversions).
//

VSubSet::VSubSet(SubSet& s) : (s)
{
  type = LINEAR_SUBSET;
  if ((holds != LINEAR_SUBSET) && (holds != NULL_SUBSET)) {
      errh.ErrorExit("VSubSet::VSubSet(SubSet&)",
                "Attempted to cast a non-vector subset to a vector subset", s);
  }
}

// ***********************************************************************
//
// Output for debugging:
//

void VSubSet::debug_out(ostream &c, int indent)
{
  static char* name = "VSubSet object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
//  Build a vector subset of a vector space.
//

VSubSet::VSubSet(char* namein, VSpace &s, GeObList &v)
{
  static char* sig = "VSubSet::VSubSet(char*, VSpace&, GeObList&)";
  Boolean errflag;

  if (v.Length() < 1) {
    errh.ErrorExit(sig, "Incorrect number of points specified", v);
  }

// Create the SubSetInfo block:

  info = new SubSetInfo;
  type = LINEAR_SUBSET;
  holds = LINEAR_SUBSET;

// Important info:

  info->t = holds;
  info->d = v.Length();
  info->s = s;
  info->accepts = s.NativeType();
  info->isOffset = -1;
  info->substart = 0;
  info->zerostart = info->d;
  info->offset = 0.0;
  this->CopyDebugName(namein);

// Make sure all the spanning objects are in the specified space:

  Matrix spantemp(v.Length(), s.Dim());
  if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
    errh.ErrorExit(sig,
           "Object cannot be mapped into specified space", v, s);
  }

// The user provides a list of geometric objects that span
// some n-dimensional portion of the m-dimensional space.  They must
// be independent.  We need to construct an orthonormal basis for
// the embedding space such that the first n vectors span the
// subset, and the remaining vectors are orthogonal to the subset.
// The orthonormal basis is built using Gram-Schmidt othogonalization.
// Normalize the first vector:

  Scalar mag = sqrt((spantemp[0] * Transpose(spantemp[0]))[0][0]);
  Matrix have = spantemp[0] / mag;

// Gram-Schmidt:

  info->tofull = this->GSO(have, spantemp, 1, info->d, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  }
  info->fullbas = this->GSO(info->tofull, IdentityMatrix(s.Dim()),
                0, s.Dim(), &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error", v, s);
  }

// Build the testing matrix and the projection matrix to convert
// "close" full space representations to subset representations:

  info->tstmtx = Inverse(info->fullbas);
  Matrix trg = ZeroMatrix(s.Dim(), info->d);
  for (int i = 0; i < info->d; i++) {
    trg[i][i] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}

// ***********************************************************************
// ***********************************************************************
//
// ASubSet Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Builds a standard affine subset for a space.  CAUTION:  Although this
// is a single-argument constructor, it is NOT intended to be used to
// cast a space into a subset!
//

ASubSet::ASubSet(Space& s)
{
  int size = s.MatrixSize();

// Create the new SubSetInfo:

  info = new SubSetInfo;
  type = AFFINE_SUBSET;
  holds = AFFINE_SUBSET;

  char* tag = this->BuildTagName("Standard affine subset for ", s);
  this->CopyDebugName(tag);
  this->DestroyTagName(tag);

  info->t = holds;
  info->d = (s.Holds() == VEC_SPACE) ? s.Dim() - 1 : s.Dim();
  info->s = s;
  info->tstmtx = IdentityMatrix(size);
  info->fullbas = IdentityMatrix(size);
  info->isOffset = size - 1;
  info->substart = 0;
  info->zerostart = size;
  info->offset = 1.0;
  info->fromfull = IdentityMatrix(size);
  info->tofull = IdentityMatrix(size);
  info->accepts = s.NativeType();
}

// ***********************************************************************
//
// Used to create a consistent standard affine subset, given a
// projective map.
//

ASubSet::ASubSet(ASubSet& a, ProjectiveMap& m)
{
  static char* sig = "ASubSet::ASubSet(ASubSet&, ProjectiveMap&)";
  Boolean errflag;

// Create the new SubSetInfo block:

  info = new SubSetInfo;
  type = AFFINE_SUBSET;
  holds = AFFINE_SUBSET;

  info->t = holds;
  info->s = m.Range().EmbeddingSpace();

  int dim = info->s.MatrixSize();

  info->d = (info->s.Holds() == VEC_SPACE) ? info->s.Dim() - 1 : info->s.Dim();
  info->accepts = info->s.NativeType();
  info->isOffset = dim - 1;
  info->substart = 0;
  info->zerostart = dim;

  char* tag = this->BuildTagName("Standard affine subset for ", info->s);
  this->CopyDebugName(tag);
  this->DestroyTagName(tag);

// The map had better be invertible.  Send the subset basis through
// the map, and recreate the subset:

  if (!m.Invertible()) {
    errh.ErrorExit(sig, "Requires invertible map", a, m);
  }

  if (a.EmbeddingSpace() != m.Domain().EmbeddingSpace()) {
    errh.ErrorExit(sig, "Domain mismatch", a, m);
  }

// This approach works, but is not elegant.
// Map the offset vector as a linear functional to keep it
// perpendicular to tangent space.  Map the subset basis
// as usual:

  Matrix pt = ((a.AugBasis())[a.IsOffset()] *
          Inverse(Transpose(m.MatrixRep())));
  Scalar mag = sqrt((pt * Transpose(pt))[0][0]);
  info->offset = a.Offset() / mag;

  Matrix span2 = a.AugBasis() * m.MatrixRep();

// Keep the offset vector intact by making it the first vector
// in the new basis:

  span2[a.IsOffset()] = span2[0];
  span2[0] = pt[0];

// Build a set of orthonormal vectors for the subspace, then
// augment with vectors orthogonal to the subspace:

  Matrix have = span2[0] / mag;
  info->tofull = this->GSO(have, span2, 1, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error A", a, m);
  }

// Swap unchanged offset vector back into position:

  Matrix hold = info->tofull[0];
  info->tofull[0] = info->tofull[info->isOffset];
  info->tofull[info->isOffset] = hold[0];

  info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
                0, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error B", a, m);
  }

// Calculate the testing matrix:

  info->tstmtx = Inverse(info->fullbas);

// Calculate fromfull:

  Matrix trg = ZeroMatrix(dim, dim);
  for (int i = 0; i < dim; i++) {
    trg[i][i] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// A subset can only be cast down to an affine subset if it is one (no
// conversions).

ASubSet::ASubSet(SubSet& s) : (s)
{
  type = AFFINE_SUBSET;
  if ((holds != AFFINE_SUBSET) && (holds != NULL_SUBSET)) {
      errh.ErrorExit("ASubSet::ASubSet(SubSet&)",
               "Attempted to cast a non-affine subset to an affine subset", s);
  }
}

// ***********************************************************************
//
// Memberwise assignment:
//

ASubSet& ASubSet::operator=(ASubSet& s)
{
  info = s.info;
  holds = s.holds;
  return *this;
}

// ***********************************************************************
//
// Output for debugging:
//

void ASubSet::debug_out(ostream &c, int indent)
{
  static char* name = "ASubSet object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
//  Build an affine subset of a vector or affine space.
//

ASubSet::ASubSet(char* namein, Space &s, GeObList &v)
{
  static char* sig = "ASubSet::ASubSet(char*, Space&, GeObList&)";
  Boolean errflag;

  if ((s.Holds() != VEC_SPACE) && (s.Holds() != AFF_SPACE)) {
    errh.ErrorExit(sig, "Vector or affine space required", s);
  }

  if (v.Length() < 1) {
    errh.ErrorExit(sig, "Incorrect number of points specified", v);
  }

// Create the new SubSetInfo block:

  info = new SubSetInfo;
  type = AFFINE_SUBSET;
  holds = AFFINE_SUBSET;

// Important info:

  info->t = holds;
  info->d = v.Length() - 1;
  info->s = s;
  info->accepts = s.NativeType();
  info->isOffset = info->d;
  info->substart = 0;
  info->zerostart = info->d + 1;
  this->CopyDebugName(namein);

// Make sure all the spanning objects are in the specified space:

  int dim = s.MatrixSize();
  int num = v.Length();

  Matrix spantemp(num, dim);
  if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
    errh.ErrorExit(sig, "Object cannot be mapped into specified space", v, s);
  }

// The subset basis will be a frame for the affine subset, so we need
// to derive vectors in the tangent subset.  Keep the last point
// in the list as the point origin for the frame, and keep a separate
// copy of the point's matrix rep. to use for finding the offset.

  Matrix span2(num, dim);
  Matrix pt = spantemp[num - 1];
  int slot = 0;
  for (int i = 0; i < num; i++) {
    if (i == num - 1) {
      span2[i] = spantemp[i];
    } else {

// The following line of code produces a bug with g++ 1.37.1
// (no problem with CFRONT 1.2) when using dynamically allocated matrices:
//    span2[slot++] = (spantemp[i] - pt)[0];
// Replace it with:

      Matrix temp = (spantemp[i] - pt[0]);
      span2[slot++] = temp[0];
    }
  }

// Build a set of orthonormal vectors for the subspace, then
// augment with vectors orthogonal to the subspace:

  Scalar mag = sqrt((span2[0] * Transpose(span2[0]))[0][0]);
  Matrix have = span2[0] / mag;
  info->tofull = this->GSO(have, span2, 1, num, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  }

  info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
                0, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error", v, s);
  }

// Testing matrix:

  info->tstmtx = Inverse(info->fullbas);

// Record the offset of the affine hyperplane (the projection of
// the original frame origin point onto the offset dimension):

  info->offset = (pt * Transpose(info->tofull[info->isOffset]))[0][0];

// Calculate fromfull:

  Matrix trg = ZeroMatrix(dim, num);
  for (int j = 0; j < num; j++) {
    trg[j][j] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}

// ***********************************************************************
//
// Build an affine subset of a projective space. The key here is
// that the combination of the "regular" point and the
// "points at infinity" specify a projective subspace;
// subtracting the points at infinity turns this into an
// affine subset.
//

ASubSet::ASubSet(char* namein, PSpace& s, GeOb& v, GeObList& inf)
{
  static char* sig = "ASubSet::ASubSet(char*, PSpace&, GeOb&, GeObList&)";
  Boolean errflag;

// The number of points at infinity that are provided needs to be
// appropriate.  For an m-dimensional affine subset, we need m+1
// projective points, of which m are designated to be points at
// infinity, and only 1 is a "regular" point.

  if (inf.Length() < 1) {
    errh.ErrorExit(sig, "Incorrect number of points specified", inf);
  }

// Create a new SubSetInfo block:

  info = new SubSetInfo;
  type = AFFINE_SUBSET;
  holds = AFFINE_SUBSET;

// Important info:

  info->t = holds;
  info->d = inf.Length();
  info->s = s;
  info->accepts = s.NativeType();
  info->isOffset = info->d;
  info->substart = 0;
  info->zerostart = info->d + 1;
  this->CopyDebugName(namein);

// Make sure all the spanning objects are in the specified space:

  int dim = s.MatrixSize();

  Matrix vtemp(1, dim);
  Matrix inftemp(inf.Length(), dim);
  GeObList vlst(v);

  if (!(this->ConvertList(vlst, info->accepts, s, &vtemp))) {
    errh.ErrorExit(sig, "Object cannot be mapped into specified space", v, s);
  }
  if (!(this->ConvertList(inf, info->accepts, s, &inftemp))) {
    errh.ErrorExit(sig,
           "Object cannot be mapped into specified space", inf, s);
  }

// Combine the matrix reps. of the two lists into a single matrix.

  int num = info->d + 1;
  Matrix span(num, dim);
  span[num - 1] = vtemp[0];
  for (int i = 0; i < num - 1; i++) {
    span[i] = inftemp[i];
  }

// Build a set of orthonormal vectors for the subset, then
// augment with vectors orthogonal to the subset:

  Scalar mag = sqrt((span[0] * Transpose(span[0]))[0][0]);
  Matrix have = span[0] / mag;
  info->tofull = this->GSO(have, span, 1, num, &errflag);
  if (errflag) {
    errh.ErrorExit(sig,
           "Objects to define subset must be independent", v, inf);
  }
  info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
                0, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error", v, inf);
  }

// Build the testing matrix:

  info->tstmtx = Inverse(info->fullbas);

// Set the offset of the affine hyperplane at 1.0 (it can be any
// arbitrary non-zero value for affine subsets of projective spaces):

  info->offset = 1.0;

// Calculate fromfull:

  Matrix trg = ZeroMatrix(dim, num);
  for (int j = 0; j < num; j++) {
    trg[j][j] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}

// ***********************************************************************
// ***********************************************************************
//
// PSubSet Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Builds a standard projective subset for a space.
//

PSubSet::PSubSet(Space& s)
{
  static char* sig = "PSubSet::PSubSet(Space&)";
  int size = s.MatrixSize();

  if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
    errh.ErrorExit(sig, "Vector or projective space required", s);
  }

// Create a new SubSetInfo block:

  info = new SubSetInfo;
  type = PROJECTIVE_SUBSET;
  holds = PROJECTIVE_SUBSET;

  char* tag = this->BuildTagName("Standard proj. subset for ", s);
  this->CopyDebugName(tag);
  this->DestroyTagName(tag);

  info->t = holds;
  info->d = (s.Holds() == VEC_SPACE) ? s.Dim() - 1 : s.Dim();
  info->s = s;
  info->tstmtx = IdentityMatrix(size);
  info->fullbas = IdentityMatrix(size);
  info->isOffset = -1;
  info->substart = 0;
  info->zerostart = size;
  info->offset = 0.0;
  info->fromfull = IdentityMatrix(size);
  info->tofull = IdentityMatrix(size);

  info->accepts = s.NativeType();
  if (info->accepts == VECTOR) {
    info->accepts = VEC_EC;
  } else if (info->accepts == AFF_VECTOR) {
    info->accepts = AFF_VEC_EC;
  }
}

// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// A subset can only be cast down to a proj. subset if it is one (no
// conversions).

PSubSet::PSubSet(SubSet& s) : (s)
{
  type = PROJECTIVE_SUBSET;
  if ((holds != PROJECTIVE_SUBSET) && (holds != NULL_SUBSET)) {
      errh.ErrorExit("PSubSet::PSubSet(SubSet&)",
       "Attempted to cast a non-projective subset to a projective subset", s);
  }
}

// ***********************************************************************
//
// Memberwise assignment:
//

PSubSet& PSubSet::operator=(PSubSet &s)
{
  info = s.info;
  holds = s.holds;
  return *this;
}

// ***********************************************************************

void PSubSet::debug_out(ostream &c, int indent)
{
  static char* name = "PSubSet object\n";
  this->data_out(c, indent, name);
  return;
}

// ***********************************************************************
//
//  Build a projective subspace, without base point specification.
//

PSubSet::PSubSet(char* namein, Space &s, GeObList &v)
{
  static char* sig = "PSubSet::PSubSet(char*, Space&, GeObList&)";
  Boolean errflag;

  if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
    errh.ErrorExit(sig, "Vector or Projective space required", s);
  }

  if (v.Length() < 1) {
    errh.ErrorExit(sig, "Incorrect number of points specified", v);
  }

// Create a new SubSetInfo block:

  info = new SubSetInfo;
  type = PROJECTIVE_SUBSET;
  holds = PROJECTIVE_SUBSET;

// Important info:

  info->t = holds;
  info->d = v.Length() - 1;
  info->s = s;
  info->isOffset = -1;
  info->substart = 0;
  info->zerostart = info->d + 1;
  info->offset = 0.0;
  this->CopyDebugName(namein);

  info->accepts = s.NativeType();
  if (info->accepts == VECTOR) {
    info->accepts = VEC_EC;
  } else if (info->accepts == AFF_VECTOR) {
    info->accepts = AFF_VEC_EC;
  }

// Make sure all the spanning objects are in the specified space:

  int dim = s.MatrixSize();
  int num = v.Length();

  Matrix spantemp(num, dim);
  if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
    errh.ErrorExit(sig,
           "Object cannot be mapped into specified space", v, s);
  }

// Build a set of orthonormal vectors for the subspace, then
// augment with vectors orthogonal to the subspace:

  Scalar mag = sqrt((spantemp[0] * Transpose(spantemp[0]))[0][0]);
  Matrix have = spantemp[0] / mag;
  info->tofull = this->GSO(have, spantemp, 1, num, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  }
  info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
                0, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error", v, s);
  }

// Build testing matrix and fromfull matrix:

  info->tstmtx = Inverse(info->fullbas);
  Matrix trg = ZeroMatrix(dim, num);
  for (int i = 0; i < num; i++) {
    trg[i][i] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}

// ***********************************************************************
//
// Build a projective subset of a projective space by selection of
// base points.  In combination, base points and other points
// define a projective subset.  The subsequent removal of the
// base points reduces the dimensionality of the subset.  The
// resulting dimensionality equals that of the corresponding
// primitive geometric form (e.g. pencil or bundle).
//

PSubSet::PSubSet(char* namein, Space& s, GeObList& bp, GeObList& v)
{
  static char* sig = "PSubSet::PSubSet(char*, Space&, GeObList&, GeObList&)";
  Boolean errflag;
  int num = v.Length() + bp.Length();

  if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
    errh.ErrorExit(sig, "Vector or Projective space required", s);
  }

  if ((v.Length() < 1) || (bp.Length() < 1)) {
    errh.ErrorExit(sig, "Incorrect number of points specified", bp, v);
  }

// Create a new SubSetInfo block:

  info = new SubSetInfo;
  type = PROJECTIVE_SUBSET;
  holds = PROJECTIVE_SUBSET;

// Important info:

  info->t = holds;
  info->d = v.Length() - 1;
  info->s = s;
  info->isOffset = -1;
  info->substart = bp.Length();
  info->zerostart = num;
  info->offset = 0.0;
  this->CopyDebugName(namein);

  info->accepts = s.NativeType();
  if (info->accepts == VECTOR) {
    info->accepts = VEC_EC;
  } else if (info->accepts == AFF_VECTOR) {
    info->accepts = AFF_VEC_EC;
  }

// Make sure all the spanning objects are in the specified space:

  int dim = s.MatrixSize();

  Matrix spantemp(v.Length(), dim);
  Matrix bptemp(bp.Length(), dim);

  if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
    errh.ErrorExit(sig,
           "Object cannot be mapped into specified space", v, s);
  }
  if (!(this->ConvertList(bp, info->accepts, s, &bptemp))) {
    errh.ErrorExit(sig,
           "Object cannot be mapped into specified space", bp, s);
  }

// Build a set of orthonormal vectors that span the space of base points.
// Then build an orthogonal basis to represent the subset, and finally
// add enough vectors to fill the basis out.  Note that we
// CANNOT build a matrix that converts basis std. coords. to full space
// std. coords (tofull), since fromfull is not invertible.

  Scalar mag = sqrt((bptemp[0] * Transpose(bptemp[0]))[0][0]);
  Matrix have = bptemp[0] / mag;
  Matrix temp = this->GSO(have, bptemp, 1, bp.Length(), &errflag);
  if (errflag) {
    errh.ErrorExit(sig,
           "Objects to define subset must be independent", s, bp, v);
  }
  temp = this->GSO(temp, spantemp, 0, num, &errflag);
  if (errflag) {
    errh.ErrorExit(sig,
           "Objects to define subset must be independent", s, bp, v);
  }
  info->fullbas = this->GSO(temp, IdentityMatrix(dim), 0, dim, &errflag);
  if (errflag) {
    errh.ErrorExit(sig, "Unexpected error", s, bp, v);
  }

// Build testing matrix and fromfull matrix:

  info->tstmtx = Inverse(info->fullbas);

  Matrix trg = ZeroMatrix(dim, v.Length());
  int slot = 0;
  for (int i = bp.Length(); i < num; i++) {
    trg[i][slot++] = 1.0;
  }
  info->fromfull = info->tstmtx * trg;
}


// ***********************************************************************
// ***********************************************************************
//
// Create the error information block:
//

static SubSetInfo errInfo("Uninitialized subset", 0);

// ***********************************************************************